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We study the dynamics in the neighborhood of simple and double unstable periodic orbits 
in a rotating 3D autonomous Hamiltonian system of galactic type. In order to visualize the 
four dimensional spaces of section we use the method of color and rotation. We investigate the 
structure of the invariant manifolds that we found in the neighborhood of simple and double 
unstable periodic orbits in the 4D spaces of section. We consider orbits in the neighborhood of 
the families xlv2, belonging to the xl tree, and the z-axis (the rotational axis of our system). 
Close to the transition points from stability to simple instability in the neighborhood of the 
bifurcated simple unstable xlv2 periodic orbits we encounter the phenomenon of stickiness as 
the asymptotic curves of the unstable manifold surround regions of the phase space occupied 
by rotational tori existing in the region. For larger energies, away from the bifurcating point, 
the consequents of the chaotic orbits form clouds of points with mixing of color in their 4D 
representations. In the case of double instability, close to xlv2 orbits, we find clouds of points in 
the four dimensional spaces of section. However, in some cases of double unstable periodic orbits 
belonging to the z-axis family we can visualize the associated unstable eigensurface. Chaotic 
orbits close to the periodic orbit remain sticky to this surface for long times (of the order of a 
Hubble time or more) . Among the orbits we studied we found those close to the double unstable 
orbits of the xlv2 family having the largest diffusion speed. The sticky chaotic orbits close to the 
bifurcation point of the simple unstable xlv2 orbit and close to the double unstable z-axis orbit 
we have examined, have comparable diffusion speeds. These speeds are much slower than of the 
orbits in the neighborhood of xlv2 simple unstable periodic orbits away from the bifurcating 
point, or of the double unstable orbits of the same family having very different eigenvalues along 
the corresponding unstable eigendirections. 



Keywords: Chaos and Dynamical Systems, 4D surfaces of section, Invariant manifolds, Stickiness 
in Chaos, Galactic Dynamics 
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1. Introduction 

The aim of this paper is to study the dynamics in the neighborhood of unstable periodic orbits in a 
rotating three dimensional (hereafter 3D) autonomous Hamiltonian system of galactic type. Instability of 
periodic orbits in such a dynamical system is of three kinds: simple, double and complex instability. The 
case of complex instability was studied in a recent paper by Katsanikas et al (Katsanikas et al. 2011). 
Here we study the two other cases. Furthermore, we examine the stickiness effects in our 3D system. This 
is the first time that the stickiness phenomenon is explicitly investigated in association with these kinds 
of instability (simple and double) in a 3D autonomous Hamiltonian system. For this purpose we consider 
the intersections of the orbits by a space of section. In a 3D autonomous Hamiltonian system the space of 
section has 4 dimensions. So, for the visualization of the 4D spaces of section it is convenient to use the 
method of color and rotation (Patsis &: Zachilas 1994). Rotation is used for understanding the geometry 
of a 3D projection and color is used to visualize the 4th dimension. 

3D rotating autonomous Hamiltonian systems are of special importance in Galactic Dynamics. They 
represent in a realistic way rotating galaxies, in which rotation plays a key role in the dynamics on the 
equatorial plane, as well as in the third dimension. The most important family for the vertical morphology 
of rotating disk galaxies is a family called xlvl by (Patsis et al. 2002). The xlvl family is introduced at 
the transition from stability to simple instability (S — > U) as the stability of the central, planar family xl 
undergoes a S — > U — > S transition within a small Jacobi integral interval at the vertical 2:1 resonance 
region (for a definition of resonances in rotating galactic potentials see Contopoulos pages 379 and 422). 
At the U — > S transition, another family of periodic orbits (hereafter p.o.), xlv2, is bifurcated from xl as 
simple unstable (Skokos et al. 2002). For larger Jacobi constants xlv2 becomes double unstable. Thus, it 
possesses both kinds of instabilities we want to study. In this paper we focus on xlv2 and we examine the 
structure of the phase space in its neighborhood, when it is simple or double unstable. However, in order 
to investigate in depth the dynamics in the neighborhood of double unstable periodic orbits we will also 
study referring to z-axis orbits. 

The importance of the dynamical phenomena at the 2:1 vertical resonance is due to the role of the 
associated families in shaping the vertical morphology of the 3D disks. They form to a large extent the 
peanut-shaped structure, which is characteristic of many disk galaxies observed edge-on (Patsis et al 2002). 
In addition, chaotic motion in the neighborhood of unstable p.o.'s has been proven to be significant for 
the dynamics of rotating galactic disks on their equatorial planes (for a review see e.g. Contopoulos 2009). 
In that respect the knowledge of the structure of phase space in the neighborhood of 3-dimensional p.o. 
in galactic systems is the first step for the investigation of a possible contribution of chaotic orbits to the 
observed vertical profiles of this type of galaxies. 

The space of section that we use in this paper for investigating the orbital behavior close to xlv2 orbits 
is defined by y = and we consider intersections with y > 0. Hereafter by "intersections" we will refer to 
intersections, which fulfill these properties, until we start discussing the dynamics close to z-axis. Then we 
will change space of section. Every consequent in our cross section is determined by 4 Cartesian coordinates, 
e.g. (x, x, z, z). In the method of color and rotation we project the consequents in a 3D subspace and observe 
them from different viewpoints. In the text the viewpoints are given in spherical coordinates. The unit for 
the distance d of the consequents of the surface of section from the observer corresponds to the longest side 
of the bounding box of the figure (more details can be found in Katsanikas and Patsis 2010). In practice 
we rotate the 3D figure in the computer screen in order to find the most appropriate viewpoint for the 
presentation on the 2D surface of the pages of a paper. We select the viewpoint that represents best the 
morphology of the structures formed by the consequents in the 3D projections of the 4D space of section 
in order to facilitate the reader to follow our description. Then we color every point according the value 
of its 4th coordinate. For this reason we use the "Mathematica" package and we map the interval of the 
values of the fourth coordinate into [0, 1]. 

• A well defined structure in all 3D projections and a smooth color variation along them is evidence for 
regular motion. 

• Mixing of colors indicates chaotic behavior of the consequents in 4D (Patsis Sz Zachilas 1994, Katsanikas 
k Patsis 2010) 
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1.1. Stability of periodic orbits 

The linear stability of the periodic orbits is studied, following the analysis by Broucke (Broucke 1969) in the 
case of the elliptic restricted three body problem and the method used by Hadjidemetriou (Hadjidemetriou 
1975) in the three body problem. The method has been also described in a recent earlier paper of us in 
the series of papers studying the dynamics in the neighborhood of p.o. in 3D rotating galactic potentials 
(Katsanikas & Patsis 2010). The characteristic equation in a 3D autonomous hamiltonian system is of the 

form A 4 + aA 3 + /3A 2 + aA + 1 = 0. and has 4 eigenvalues Ai, 1/A 1; A2, I/A2, where \, j- = bl± V b i 

i = 1,2, with 6i j2 = a± ^& anc [ A = a 2 - 4(/3 - 2). Thus, the stability of the p.o. depends on the two 
stability indices b\ and 62 and on the quantity A. 

For A > 0, \b\\ < 2 and |6 2 | < 2, all four eigenvalues are on the unit circle in the complex plane and 
the orbit is called stable (S) (see e.g. Contopoulos 2002, Contopoulos & Magnenat 1985). 

Then we have three kinds of instability: 

(1) simple instability (U): A > and > 2, I&2I < 2 or |&i| < 2, |6 2 | > 2, (two eigenvalues on the unit 
circle and two on the real axis) 

(2) double instability (DU): A > and |&i| > 2 and I62I > 2 (four eigenvalues on the real axis) 

(3) complex instability (A) (all four complex eigenvalues are off the unit circle in the complex plane): 
A < 

(for details see Contopoulos 2002, p. 286). Hereafter we will indicate the periodic orbits as S, U, DU and 
A according to their stability. A generalization of the three kinds of instability in Hamiltonian systems of 
N > 3 degrees of freedom can be found in Skokos (Skokos 2001). 

In section 2 we present our Hamiltonian system. In section 3 we discuss the stability of a family of p.o. 
we call xlv2. In section 4 we study the orbital behavior close to simple unstable orbits of xlv2, in section 
5 we investigate the dynamical behavior in the neighborhood of double unstable orbits of xlv2, as well as 
of the z-axis families, in section 6 we compare the diffusion of the chaotic orbits away from the p.o. and 
finally in section 7 we discuss and enumerate our conclusions. 

2. The Hamiltonian System 

The Hamiltonian of our system in Cartesian coordinates is: 

H(x, y, z, x, y, z) = ^(± 2 + y 2 + z 2 ) + $(x, y, z) - Ul 2 {x 2 + y 2 ), (1) 
where $(x, y, z) is the potential: 



* (x y Z ) = 9Mi 9Mi m 

(x 2 + $ + [a 1 + (£ + 6?)V2]2)l/2 {X 2 + £ + [o2 + + ^)l/2 ]2) l/2 

This potential in its axisymmetric form (q a = l,q b = 1) can be considered as an approximation of the 
potential of the Milky Way (Miyamoto & Nagai 1975). 

In our units, the distance R=l corresponds to 1 kpc. The velocity unit corresponds to 209.64 km/sec. 
The Jacobi constant Ej = 1 corresponds to 43950 (km/ sec) 2 and we will call it "energy". For the pa- 
rameters we have used the following values: a\ = kpc, b\ = 0.495 kpc, M\ = 2.05 x 10 10 M G , a>2 = 
7.258 kpc, b 2 = 0.520 kpc, M 2 = 25.47 x 10 10 M , q a = 1.2, q b = 0.9. For the angular velocity we used 
the value f2& = 0.045 in our units. The parameters q a , qb determine the geometry of the disks, while a, b 
are scaling factors. The specific values of the geometrical parameters do not have a particular physical 
meaning. The dynamical phenomena we study in our paper are typical for rotating triaxial systems. The 
0;, value we use here is different from the one we used in (Katsanikas & Patsis 2010) in order to have larger 
regions in the parameters' space with non-escape orbits close to double unstable periodic orbits (Section 
5). 
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3. Stability Diagram 

In order to study the evolution of the stability of a family of periodic orbits as a parameter of our system 
varies we use the stability diagram (Contopoulos 2002 p. 287). The stability diagram gives the values of the 
coefficients b\ and 62 for successive periodic orbits in a family by changing only one parameter of the system 
(in our case Ej), while keeping the other parameters constant. In Fig. Q] we observe in such a diagram how 
the family of periodic orbits xlv2 (Skokos et al. 2002) bifurcates from the 2D central family xl on the plane 
z = at the U — > S transition that happens in the vertical 2:1 resonance region (Contopoulos 2002 p. 390). 
This bifurcation occurs for Ej ~ —4.7. We observe that the stability parameters b\, 62 of the family xlv2 
decrease as Ej increases. The value of 62 is initially about 1.95 at Ej = —4.7 (left side of diagram). When 
62 becomes —2 we have a transition of xlv2 from simple to double instability (U — > DU) (for Ej ~ —3.41) 
and the xlv2 periodic orbits become double unstable. The Ej interval of xlv2, where it is simple unstable 
(—4.7 < Ej < —3.41) and the Ej interval of xlv2 where it is double unstable (—3.41 < Ej < —2.3), will 
be called hereafter regions U and DU respectively (Fig. [p. 

4. Simple Instability 

First we study the orbital behavior of the xlv2 family in the region U. In the following applications we 
use the set of projections (x,x,z) and (x,z, z). However, we find qualitatively similar results by using the 
(x, x, z) and (x, z, z) projections respectively. The use of the two projections gives us all needed information 
to understand the structure of the distribution of the consequents in the 4D space by means of the color- 
rotation method, as we will see below. 

The initial conditions to be studied are chosen in two different ways. First we investigate the dynamical 
behavior of the orbits in the neighborhood of the simple unstable periodic orbit xlv2 at Ej = —4.66 by 
varying the values and the direction (Ax, Ax, Az or Az) of the perturbation. Secondly, we consider 
a constant perturbation Ax = 10 -4 and we apply this perturbation to the initial conditions of several 
periodic orbits of the xlv2 family in the region U for different Ej values. We start by keeping Ej constant 
and investigate as a sample case the dynamical behavior very close to the point where we have the U — > S 
transition of the xl family and the bifurcation of the family xlv2 as U. 

4.1. Orbital behavior close to the bifurcating point - A sample case 

The energy Ej = —4.66 is indicated with an arrow in Fig. [U The periodic orbit has initial conditions 
(xq, xq, zq, zq) = (0.037941612, 0, 0, 0.42902353). In order to study the orbital behavior in its neighborhood 
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we added a small perturbation in the i-direction equal to Ai = —6 x 10 . We call this perturbed orbit 
"u". The first 2580 consequents of orbit "u" define a ribbon of the form of a double loop in the (x,x,z) 
subspace (Fig. [2]). In the bottom of the figure we give also from left to right the 2D projections (x, x), (x, z) 
and (x, z) in order to help the reader understand the topology of the structure. The two lobes of the double 
loop are inclined with respect to each other and form an angle around 50° in the (x, z) projection (middle 
panel in the lower part of Fig. [2]). The (x, z) projection clearly also shows that the thickness of the lobes is 
very small (but not zero) and this is the reason we describe them as ribbons. Essentially the branches of the 
lobes are thin tubes. The 4D representation of this structure, colored according to the normalized z values 
and using (x, x, z) for the 3D spatial projection is given in Fig. [2] (top). For the 3D projection we use here 
a point of view defined by the angles (9, <fi) = (30°, 108°). We observe that this ribbon has a smooth color 
variation from red to violet, in agreement with structures found in the neighborhood of simple periodic (U) 
orbits by Patsis and Zachilas (1994). The structure of the phase space close to simple unstable periodic 
orbits has been studied for the first time by Magnenat (1982). However, Magnenat used 2D projections and 
these did not allow the full understanding of the 4D structure of the double loop objects. In the present 
case, our 4D representation allows us to understand the details of this structure. 

By inspection of Fig. [2] we infer that the color variation along the two branches is such that at the 
intersection region we have the meeting of two different shades. This means that in this region the points 
have different values in the 4th dimension and there is no real intersection in the 4D space of section. 
Consequently the ribbon-like double loop structure has not an intersection in 4D. We observe that the 
colors at the region of the intersection are not only different, but that they are at the two extremes of our 
color scale, i.e. violet and red (see right side of Fig. [2]). The colors in this region are clustered around two 
values (i = ±0.429) very close to the z initial condition of the periodic orbit z = 0.42902353. 

Figure [2] reflects the symmetry of the galactic potential we study. Due to this symmetry, the families 
of p.o. that bifurcate at the vertical resonances from the planar xl family come in pairs with mirror 
morphology with respect to the equatorial plane of the system in the configuration space, i.e. xlvl and 
xlvl', xlv2 and xlv2' etc. The p.o. of xlv2 and xlv2' at every Ej differ only in the sign of their zq initial 
conditions. In the 3D (x, x, z) projection we use in Fig. [2j the initial conditions of xlv2 and xlv2' are 
identical and the p.o. is in the center of the intersection region of the two branches of the double loop 
ribbon (indicated with an arrow). The clustering of the colors around the two values (zq ~ ±0.429) in the 
central region is due to the coexistence of xlv2 and xlv2'. The consequents drift away from the central 
region of the double loop by deviating from the two (±io) initial conditions with different colors. The first 
2580 consequents are along the double loop structure, thus they indicate a regular rather than chaotic 
motion. 

We reach the same conclusion by using other 3D projections for the spatial coordinates of the con- 
sequents. In Fig. [3] the first 2580 consequents of the orbit "u" in the 3D projection (x,z, z) form again a 
ribbon-like surface, however this time with a single loop morphology. Color is determined by the x value 
of the consequents and we observe a smooth color variation along the ribbon. In this case the two repre- 
sentatives of the p.o xlv2 and xlv2' are at two different positions in the (x, z, z) space, since they have 
different zq initial conditions. Their positions are indicated with arrows in Fig. [3l The regions close to both 
of them are colored red (same x values) and departing from them, red becomes orange, then yellow, green, 
light blue, blue and finally violet. The ribbon is a 4D object. From both Figs. [2] and [3] we conclude that 
the first 2580 consequents form a well-defined ribbon-like surface that has no self-intersections in 4D. 

A similar analysis can be done by means of the two other 3D projections (x, x, z) and (x, z, z), which 
give 4D structures morphologically similar to (x, x, z) and (x, z, z) respectively. It is evident that the orbit 
"u" in the neighborhood of a U xlv2 p.o. for 2580 consequents is confined in a particular region of the 
space of section with dimensions (Ax x Ax x Az x Ai) < (0.00035 x 0.0005 x 0.04 x 1). As we will see 
below, for larger integration times the orbit diffuses to a larger region of the phase space. The structures 
we find in the phase space indicate a sticky behavior and this will be further analyzed in section 4.4. 

We note that although in general we have different colors in the intersections of double loop type 
structures, like the one in Fig. [2h., in exceptional cases we can have also real intersections in the 4D space. 
Such a case is encountered in the neighborhood of U z-axis orbits in 4D spaces of section defined by 
z = 0. The z-axis orbits lie entirely on the rotational, z axis, of the system. They have initial conditions 
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Fig. 2. Upper part: A 4D representation of the orbit "u" for 2580 consequents. The consequents are given in the (x,x,z) 
projection, while every point is colored according its z value. We observe that an double loop ribbon is formed with a smooth 
color succession on it (cf. with the color bar on the right). The location of the periodic orbit (p.o.) is indicated by an arrow. 
Our point of view in spherical coordinates is (9,cf>) = (30°, 108°). Lower part: The three 2D projections (x, x,), (x, z) and 
(x, z). The ranges on the axes are as in the main figure. 



xlv2 




Fig. 3. 2580 consequents of the orbit "u" in the 3D subspace (x,z,z). This time every point is colored according its x 
value. The arrows indicates the locations of the periodic orbits xlv2 and xlv2'. Our point of view in spherical coordinates is 
(e,<j>) = (45°, 144°). 

(x, x, y, y) = (0, 0, 0, 0) and are discussed in detail in section 5.2 in the context of the study of double 
instability. 

If we compute a larger number of intersections we observe a drastic change in the topology of the 
ribbon-like surface. The points deviate from the well-defined 4D structures that we described above and 
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Fig. 4. Mixing of colors for the first 1000 consequents of the orbit "u" (beyond the first 2580 consequents) that deviate from 
the 4D ribbon surface. The 3580 points totally are depicted in the 3D subspace (x, z, z) and the color represents the 4th 
dimension x of the points. Our point of view in spherical coordinates is (0, <j>) = (45°, 144°). 

form a cloud of points around them in the 3D projections with mixed colors in their 4D representations. 
However, in the particular case we study here, for about 3580 intersections, i.e. for 1000 consequents 
more than those building the structures in Figs. [2] and El in the 3D projections that include the z and z 
coordinates we discern a ® structure (i.e. a greek letter O rotated by 90°) surrounded by a ring. This can 
be seen in Fig. EJ where we depict the consequents of the orbit "u" in a (x, z, z) projection and color them 
according to their x values. This new structure in the 3D projection of the phase space is maintained for 
several thousands of intersections more. We have again a confinement of the consequents in a restricted 
volume of the phase space, which is larger than previously. After 3580 intersections the consequents extend 
within a region (Ax x Ax x Az x Ai) < (0.18 x 0.7 x 0.8 x 2.4). In the new situation we have on one hand 
the appearance of structures in certain projections which indicates order, and on the other hand mixing of 
colors, which indicates chaoticity. 

The origin of the observed structure can be understood best in the (z, z) projection (Fig. [5]). We remind 
that we study the structure of phase space in the region U, which starts with the bifurcation of the family 
xlv2 at the U —> S transition of xl in the vertical 2:1 resonance region (Fig. [1]). For smaller energies in 
the resonance region, at Ej = —5.028, we have the S —> U transition of xl and the introduction in the 
system of the stable family xlvl (and its symmetric family xlvl'). Thus for Ej = —4.66, we have, apart 
from the two simple unstable p.o. xlv2 and xlv2', the two stable xlvl and xlvl' and the planar stable 
xl. The presence of all these simple periodic orbits determines the orbital behavior of "u" in the following 
way: 

In the central part of Fig. [5h. we observe a red "ellipse with corners" , which forms within the first 2580 
intersections and which is the projection of the structure of Fig; [3] in the (z, z) plane. The initial conditions 
of xlv2 and xlv2' are located at the upper and lower corner of the red elliptical curve respectively. The 
blue ellipse inside the red one is the projection of a rotational torus (Vrahatis et al. 1997) belonging to a 
quasiperiodic orbit found by perturbing xl in the z direction. During the first 2580 consequents the orbit 
"u" stays sticky to the outer 4D rotational tori of xl. The outermost rotational tori belonging to xl are 
broader than the ribbon-like tori of the "u" orbit. In Fig. [6] we give two characteristic projections from 
the point of view (0,(f>) = (90°, 42°) in spherical coordinates. In (a) we have the (x,x,z) and in (b) the 
(x, z, z) projection. In both of them we overplot the 2580 consequents of "u" (red) and one of the outermost 
rotational tori around xl (black). 

By continuing the integration of the orbit "u" for more than 2580 intersections we see that the next 
550 intersections form a CD-like structure with two lobes on the sides of the red elliptical structure (Fig. [5^, 
indicated with a couple of white arrows). These consequents surround the rotational tori around the stable 
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Fig. 5. (a) 5000 consequents of the orbit "u" in the 2D (z,z) projection. The first 2580 consequents form the red "ellipse 
with corners" in the central part of the figure. The location of xlv2 and xlv2' are indicated with black symbols on it (up and 
down). The blue ellipse inside the red one is the projection of a rotational torus around xl. The rest of the consequents, after 
the first 2580, diffuse in phase space, but remain in the region around the rotational tori belonging to the stable p.o. of the 
xlvl and xlvl' families. The projections of the two tori around each of the two stable p.o. (indicated also with black symbols 
in the middle of the "empty" regions) appear in this case like green and magenta invariant curves. Arrows point to the lobes 
and ring structure (see text), (b) 15000 consequents of the orbit "u" in the 2D (z,z) projection. They visit all available regions 
of the phase space. However the sticky regions formed during the integration time are clearly discernible. 



p.o. xlvl and xlvl', which are very thin and appear in this projection like invariant curves (green and 
magenta). The rotational tori of these two families are still in the phase space region close to the U p.o. 
xlv2, since they have been bifurcated at a nearby Ej from xl (Fig. [T]). The initial conditions of xlvl and 
xlvl' are in the middle of the rotational tori, marked with black symbols. 

Considering more than 3130 consequents we find that the next 1370 intersections stay trapped in a 
ring in the (z, z) projection, which in its turn surrounds all structures we have found in the phase space up 
to now. We point to this ring with two black arrows in Fig. [5h. (right above and left below). Between the 
two lobes surrounding the rotational tori around xlvl and xlvl' and the ring structure, we observe a chain 
of eight stability islands (colored black), which are the (z,z) projections of the corresponding rotational 
tori of a stable 8-periodic orbit. Finally, after 4500 intersections, the consequents start diffusing into a 
larger volume of the phase space. Figure [5^ depicts the space of section after 5000 intersections. After 
15000 intersections the consequents have occupied all available phase space (Fig. (5b)- However, the denser 
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Fig. 6. (a) The (x,x,z) and (b) the (x,z,z) projections of the first 2580 consequents belonging to "u"(red) and one of the 
outermost rotational tori around xl (black). The point of view is given by the angles (9,(f>) = (90°, 42°). 



regions clearly show us where the orbits have been trapped for longer times. 

We have to compare the structures described up to now with the 4D space of section we presented in 
Fig. |U The projections in Fig. 0are part of Fig. [U We observe them as we rotate the (x, z, z) diagram in 
our computer screen. Since Fig. H] consists of 3580 consequents, apart from the ribbon it includes also the 
consequents forming the CD- like structure around the rotational tori of xlvl and xlvl' as well as part of the 
consequents that form the "ring" in Fig. In Fig. [J] we emphasize on the distribution of the consequents 
in the 4th dimension. We observe that in Fig. 0] the color, representing the normalized ^-coordinate values, 
is dominated by green shades. This means that although we do not have a smooth color variation as in 
Figs. I2|3I the consequents stay confined in a restricted volume of phase space. This happens indeed. The 
x values for the first 3130 intersections are in the range —0.008 < x < 0.08, while within the next 450 
intersections this range is almost doubled. Due to this expansion of the range of x values, in Fig. [4] are 
introduced the orange-red and blue-violet consequents. The fact that in 3D projections we find structure 
and that the color variation is dominated by a particular shade, indicates a weakly chaotic behavior of 
the orbit "u" . This information is included in the 4D representation of the phase space as it is given in 
Fig. HI By combining Figs. H] and [5] on a computer screen, we reveal not only the weakly chaotic character 
of the orbit, but also the regions where it sticks during the integration time. This shows the efficiency of 
the color-rotation method. 

One can calculate for the orbit "u" the finite-time Lyapunov Characteristic Number, LCN(t), which 
is defined by 

LCN(t) = - In 

where £(io) an d are the distances between two points of two nearby orbits at times t = and t 
respectively (see e.g. Skokos 2010). As we can see in Fig. [7J the changes we have found in the orbital 
behavior of "u" do correspond to characteristic changes of the slope of the LCN{t) curve. During the first 
2580 intersections, indicated by the left arrow in Fig. [TJi the curve decreases towards zero. Then during 
the next 450 intersections, until the time pointed by the right arrow in Fig. [7^ when the consequents form 
the CD-like structure with the lobes around the xlvl and xlvl' tori, it continues to decrease but the curve 
has a much smaller slope. After that it starts increasing and its evolution after this point until i = 3x 10 4 
can be followed in Fig. [7b (note the different scales in the vertical axes). LCN{t) increases, then levels off, 
and has a decreasing part until the left arrow in Fig. [7b. which corresponds to about 4180 intersections. 
At this point we have an expansion of the consequents in phase space, which can be understood in the 
(x,x) projection (not shown in a figure). Beyond the left arrow in Fig. [7b the LCN(t) curve increases. 
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Fig. 7. The variation of LCN(t) for the orbit "u". The "ribbon" 4D morphology lasts until the point indicated with the left 
arrow in (a). The thin layer around the xlvl and xlvl' rotational tori that corresponds to the lobes (Fig. [5^) is formed within 
the time interval between the two arrows in (a). The left arrow in (b) indicates a time at which we have an expansion of the 
consequents observed in the (a:, x) projection, while to the right of the right arrow in (b) the consequents start visit all the 
available phase space, (c) Finally the curve levels off at a value LCN(t) ~ 2.74 x 10~ 2 . 



Then, between Fig. [7b and [7J; the curve decreases and in Fig. 03 it increases again. Finally, for longer 
times of integration, when the orbit explores all the available phase space, the curve levels off at a value 
LCN(t) m 2.74 x 1(T 2 (Fig. [7b) 

The right arrow in Fig. [7b corresponds to 4500 intersections and the region to the right of this arrow 
corresponds to the expansion of the consequents from the ring to the whole available volume of phase space 
as we observe in Fig. [5b. Although LCN(t) gives globally the variation of the volumes of the phase-space 
filled by the consequents, by means of the color-rotation representation we can see in detail the specific 
regions where the consequents are trapped and the coordinate along which they diffuse to occupy larger 
volumes of the phase space. This information is essential in Galactic Dynamics for building self-consistent 
models using the orbital theory. 

Summarizing the investigation of the orbital behavior of "u" in the neighborhood of a simple unstable 
periodic orbit we underline that it is characterized by two sticky periods, one associated with the rotational 
tori of xl and a second sticky period associated with the rotational tori of xlvl and xlvl'. During the 
latter sticky period we observe the ® structure and the ring in the (z, z) projection. In the following section 
we examine how general is this behavior. 

4.2. The general orbital behavior in region "U". 

A dynamical behavior similar to that of the orbit "u" is observed by perturbing the in initial condition of 
xlv2 up to Az ~ — 10~ 2 (i.e. an absolutely largest \Az\ = 10 -2 ). We find a similar behavior if we perturb 
x by at most Ax ~ 2 x 10~ 4 . Initial conditions that have the same xq, xq and zq with xlv2 and zq that 
deviates by Ai < — 10~ 2 , give rotational tori around xl (see Katsanikas & Patsis 2010). 

Another class of orbits in the neighborhood of xlv2 is characterized by immediate stickiness around 
the xlvl and xlvl' rotational tori, without building initially a "ribbon" . We find them by keeping the same 
energy constant (Ej = —4.66) and computing orbits by adding the perturbations given in Table 1 to the 
initial conditions of xlv2. In all these cases we do not observe the well defined ribbon-like surface that we 
have seen in Figs. [2] and El The orbits, starting close to the p.o. xlv2, depart directly from its immediate 
neighborhood and stick around the rotational tori of xlvl and xlvl', before they expand and explore all 
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Table 1. The range of deviations from the ini- 
tial conditions of the periodic orbit xlv2 for 
Ej = —4.66 that give orbits whose conse- 
quents stay sticky to xlvl and xlvl' and then 
diffuse occupying all the available phase space. 





positive direction 


negative direction 


Ax 


Aa; < 2 x 10 _1 


Ax > -2 x KT 1 


Ai 


2 x 10~ 4 < Ai < 1.2 


Ai > -2 


Az 


Az<7x 1CT 2 


Az > -7 x icr 2 


Az 


Az < 1.2 





available phase space. They form CD-shaped structures surrounded by rings as in Fig. [4] without the central 
elliptical structure we observe in Fig. In Table 1 we give the ranges of the perturbations of the initial 
conditions for which we have observed this orbital behavior. For deviations from the p.o. in the z direction 
in the interval 4.4 x 1CP 1 < Ai < 4.7 x 1CP 1 we reach the eight tori of the stability islands (Fig. [5]). By 
perturbing an initial condition we keep all others equal to the corresponding initial condition of the p.o. 

For perturbations larger than those in Table 1 in the x, x or i-direction we do not find any stickiness 
but the consequents form clouds around the initial conditions of the p.o. xlv2. These clouds are not 
confined in a volume smaller than the total available phase space. They occupy always larger volumes as 
the integration time increases. We have to take into account that for larger Az perturbations than the 
listed in Table 1 we first fall on initial conditions of the rotational tori of xlvl and xlvl' and then for even 
larger perturbations we find again the chaotic sea. 

The second way we followed to explore the structure of phase space around xlv2, is by varying Ej, 
keeping the same perturbation in its initial conditions. We did so by adding a Ax = 10 -4 perturbation in 
its xq initial conditions, for Ej's in the region U. For —4.7 < Ej < —4.4 the orbital dynamics of the orbits 
are as in the case for Ej = —4.66 that we described in section 4.1. However, for Ej > —4.4 a cloud of points 
surrounds the initial conditions of xlv2 in the 4D phase space. This means that in all 3D projections we 
observe clouds surrounding xlv2 from the beginning of the calculation and also a mixing of colors. 



4.3. Asymptotic Curves 

Let us now try to understand the orbital behavior of the orbits we have studied in the previous section 
by investigating the structure of the asymptotic curves of xlv2 in the 4D spaces of section. The asymp- 
totic curves are computed by mapping a line segment in the space of section on an associated dilating 
eigenvector, starting very close to a simple unstable periodic orbit xlv2 (see e.g. Magnenat 1981). The 
example we present is again for Ej = —4.66. We have already studied the structure of the phase space in 
the neighborhood of this simple unstable periodic orbit in sections 4.1 and 4.2. The initial conditions of 
the p.o. have been given in section 4.1. 

We consider the unstable eigenvector 
(xi,xi,zi,ii) = (-0.0011989,-0.00109617,-0.0459267,0.998943) and the corresponding eigenvalue A = 
1.63536. We compute the asymptotic curves of the unstable manifold and we study the morphology of its 
asymptotic curves. The initial segment is composed of the initial conditions {x 2 , x 2 , z 2 , ^2) of the form: 

x 2 = x + CX\ 
X 2 = X + cx\ 

z 2 = Z + CZ\ 

Z 2 = Z + CZi 

First we compute the asymptotic curves for c < and we call this part of the unstable manifold, part 
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Fig. 8. The first 50 consequents of a single initial condition with c = — 10~ 3 in part A of the unstable manifold in the (x, x, z) 
projection. The periodic orbit (p.o.) is indicated with a black symbol x. The red arrows point to the 3rd, 16th, 24th and 
the 38th consequent. The successive consequents (red crosses) follow the directions of the black arrows. Our point of view in 
spherical coordinates is (9,<f>) = (78°, 84°). 



A. The initial segment for part A is composed of 10 4 points (for — 1CP 3 < c < — 1CP 7 with a step 10 -7 ). For 
each point, 50 consequents were calculated. The basic criterion for choosing the values of the perturbations 
and the number of the consequents is to stay on the asymptotic curves with sufficient accuracy. The 
monodromy matrix conserves its symplectic character (determinant = 1) with an accuracy at least of the 
order of 1CT 6 . 

In Fig. [8]we observe the first 50 consequents of part A of the unstable manifold in the (x, x, z) subspace 
of the 4D space of section. The first 2 consequents are very close to the periodic orbit. Then, we discern the 
3rd consequent (indicated with an arrow). The next 13 consequents follow the direction of the black arrows, 
starting to the right and upwards, and they form a loop. After that, we return back to the neighborhood 
of the periodic orbit. Then, the 16th consequent as well as the following 7 remain again very close to 
the periodic orbit. From the 24th consequent and for the next 14 consequents we observe the formation 
of another loop, below the first one, following the directions that are indicated by the arrows. The final 
12 consequents are again very close to the periodic orbit. The evolution in time of an initial segment of 
initial conditions leads to the formation of a ribbon-like double loop structure, similar to the one in Fig. 
[2j The same holds for the (x, z, z) 3D projection of the manifold, which has a similar structure with the 
corresponding projection of the orbit "u" (Fig. [3]). As in Figs. 12131 color on the asymptotic curves varies 
smoothly, so the 4D morphology is as for the orbit "u". 

For c > in the initial segment of the initial conditions the asymptotic curves evolve in the way we 
describe below and form the second part of the unstable manifold. We call this second part, part B. The 
initial segment of initial conditions for part B is composed by 2 x 10 points (for c = 10~ 6 < c < 2 x 10 -2 
with a step 10 -6 ). For every point we computed this time 20 consequents. In Fig. [9] we depict the first 
60 points of the initial segment in the (x, x, z) projection and color them according to their z values. The 
figure helps us understand that the asymptotic curves of part B of the unstable manifold are 4D curves, 
which wind around the point corresponding to the xlv2 p.o. in the space of section. Note that we have 
a smooth color variation along the curve from light blue to red, following the direction of the successive 
numbered arrows. 

Evolving the orbits for a larger number of points (consequents) from the initial segment (400 instead 
of 60 points) we observe that the asymptotic curves of part B oscillate rather far in the (x, x, z) projection 
(Fig. [Toll . For example in the z-direction these oscillations extend between —0.4 < z < 0.4. Considering 
the evolution of all 2 x 10 4 consequents of the initial segment the asymptotic curves fill roughly the 4D 
space (x 1 ,x 2 ) x {x 1 ,x 2 ) x (z 1 ,z 2 ) x (ii,i 2 ) = (0.025,0.05) x (-0.08,0.08) x (-0.4,0.4) x (-1,1) (Fig.flOj). 
The figure is complicated and one can discern apart from the formation of loops at the upper left corner 
in Fig. [I0l a ring/disk structure that has at its center the initial conditions of the p.o. 
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Fig. 9. The 4D space of section for the evolution of the first 60 points of the initial segment of part B in the case of the 
unstable manifold. The points are depicted in the (x, x, z) subspace and are colored according their z value. The successively 
numbered arrows indicate the direction of evolution of the consequents along the manifold. The location of the periodic orbit 
(p.o.) is indicated by an arrow. Our point of view in spherical coordinates is (8, <j>) = (56°, 22°). 




Fig. 10. Part B of the unstable manifold in the 3D subspace (x, x, z) of the space of section is given by following the evolution 
of 2 x 10 4 points in the initial segment for 20 consequents (a total of 4 x 10 5 consequents). The periodic orbit is indicated by 
a black "+". Our point of view in spherical coordinates is in both panels (6, (j>) — (56°, 80°). 

As in the case of the orbits we have discussed in subsections 4.1 and 4.2, we also study the asymptotic 
curves by combining several 3D projections in which we color their points according to the fourth coordinate. 
Using the (x, z, z) 3D projection and coloring the points according to their x values the asymptotic curves 
are depicted in the following way. 

In Fig. [11] we give the first 400 consequents totally (20 points on the initial segment evolved for 20 
intersections). They form an asymptotic curve that oscillates in the 3D (x,z,z) projection. Along the 
asymptotic curve the color variation is smooth (from light blue to red to dark blue and then back to light 
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Fig. 11. The first 20 points of the initial segment of part B of the unstable manifold evolved for 20 intersections (a total of 
400 consequents) in the 4D space of section. The points are depicted in the (x, z, z) subspace and are colored according its x 
value. The arrow labeled with "p.o." indicates the location of the periodic orbit (p.o.). The consequents form a loop following 
the journey that is indicated with black arrows. Our point of view in spherical coordinates is (8, cj>) = (30°, 45°). 

blue), which means that the asymptotic curve is a 4D curve. Evolving further the asymptotic curve in time 
we observe that the oscillations become more complicated but the smooth color variation continues along 
the curve. This clearly shows that the oscillations, which have been calculated by Magnenat (Magnenat 
1982) in 2D projections, occur in the 4D space. 

In order to study the large scale evolution of part B of the unstable manifold we evolve next 400 points 
of the initial segment. The result is depicted in Fig. [12l After one oscillation on the left side of Fig. [12a 
we have several oscillations in the right part of the figure. The oscillations become longer than those of 
the previous case (Fig. [TT]) and they extend to larger absolute values of z in the 3D subspace (x, z, z) (Fig. 
[12a). By evolving 20000 points of the initial segment, the asymptotic curves form a G^shaped surface in 
the (x,z, z) projection (Fig. 112b) consisting of many parallel filaments. The filamentary character of this 
structure is due to the parallel oscillations of the asymptotic curves in the (x, z, z) subspace. 

There is an obvious similarity of Fig. [T2]with Figs. H]and[I3 The projection of the asymptotic curves 
in five of the six 2D projections show clouds of points that seem to fill densely the available phase space. 
However, in the (z, z) projection we get a clear insight of how the curves evolve in time. In this projection 
(Fig. [13}) we see that the asymptotic curves of part B of the unstable manifold depart from the p.o. in the 
unstable direction (arrow 1) and then they follow the arcs as indicated by the arrows 2,3,4,5 and 6. The 
asymptotic curves oscillate along a direction vertical to the stable eigendirection. 

The morphology we encounter for the case of the stable manifold is similar to that of the asymptotic 
curves of the unstable invariant manifold. Below we discuss the morphology of the asymptotic curves for 
other values of the Jacobi constant in the region U of simple unstable xlv2 periodic orbits. 

4.3.1. Energy variation and asymptotic curves 

We describe now the invariant manifolds in the neighborhood of a simple unstable periodic orbit for values 
of Ej —4.7 < Ej < —4.62. They have again two parts. 

As the energy increases, for values larger than Ej = —4.66, the asymptotic curves of part A form a 
ribbon, double loop, surface in the 3D subspace (x, x, z). As an example we give the case with Ej = —4.63 
(Fig. [14"|) . Although the lobes of the double loop are distorted, they form an acute angle in their (x,z) 
projection like the corresponding lobes of the orbit "u" (lower middle panel in Fig. [2]). We observe also 
that the consequents wind around the periodic orbit (as the consequents in Fig. [9|) when they approach 
the periodic orbit. The manifold has a similar structure like the one we have encountered in the central 
region of "u" in the same 3D projection. Also, in these higher energies, the asymptotic curves of part A 
form again a circular ribbon-like surface in the 3D subspace (x, z,z) similar to the orbit in Fig. [3l The 
increase of energy causes oscillations of the asymptotic curves of part A, in the 3D subspace (x, z, z). 
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Fig. 12. The asymptotic curves in the 3D subspace (x,z,i). (a) the first 400 points of the initial segment of part B of the 
unstable manifold for 20 consequents (total 8000 consequents). The periodic orbit is given with a black x symbol. After 
performing a loop to the left side of the figure the asymptotic curves make several oscillations in the central and the right part 
of the figure. The z values of the curve extend roughly within a Az — 2 range. Our point of view in spherical coordinates is 
(9, <jj) = (68°, 85°). (b) for 2 x 10 points of the initial segment of part B of the unstable manifold for 20 consequents (a total 
of 4 x 10 5 consequents). The location of the periodic orbit (p.o.) is indicated with an arrow. Our point of view in spherical 
coordinates is (0,0) = (66°, 84°). 



In the interval —4.61 < Ej < —3.47 the invariant manifolds do not have two different forms. Along 
both directions the asymptotic curves have the same morphology. The asymptotic curves resemble the 
asymptotic curves of part B of the unstable invariant manifolds for Ej = —4.66 (the case that we described 
in the section 4.3). 

4.4. Stickiness 

From the investigation of the spaces of section and the numerical calculation of the asymptotic curves 
we presented in the two previous subsections it becomes evident that the phenomenon of stickiness not 
only is present but characterizes the dynamical behavior of many orbits that start in the neighborhood 
of the simple unstable orbit xlv2. As we know, stickiness in chaos refers to chaotic orbits that stay in a 
particular region of the phase space, before escaping to large distances. In 2D autonomous Hamiltonian 
systems stickiness appears near the borders of stability islands, as well as near the asymptotic curves of 
the unstable manifolds for a long time (Contopoulos & Harsoula 2008, Contopoulos & Harsoula 2010). In 
our paper we explore this phenomenon in the 4D space of section of our 3D system and we try to explain it 
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Fig. 13. The 2D projection (z, i) of part B of the unstable manifold in the space of section for 20000 points as in Fig. 112b . 
The arrows show the direction followed by the consequents. The stable and unstable eigendirections are depicted with magenta 
and a green lines respectively. The periodic orbit xlv2 is indicated with a black "+" symbol. 
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Fig. 14. The (x,x,z) 3D projection of part A of the unstable manifold for Ej = —4.63. The asymptotic curves are computed 
for 10 consequents and for 450 initial conditions (for — 10~ 6 < c < —4.5 x 10 -4 with a step — 10~ 6 - a total of 4500 
consequents). The location of the periodic orbit is indicated with a black "x" symbol. Our point of view in spherical coordinates 
is (6,<f>) = (64°, 178°). 



in combination with the structure of the associated manifolds. It is the first time that stickiness in chaos is 
studied in detail in the neighborhood of simple unstable periodic orbits of a 3D autonomous Hamiltonian 
system. In order to study the phenomenon of stickiness in our system we consider first the representative 
example of a periodic orbit in the region U of the xlv2 family for Ej = —4.66. 

If we add a negative perturbation in the i direction (as we mentioned in section 4.2), or add a small 
positive perturbation in the x direction (for 10~ 6 < Ax < 2 x 10~ 4 ) to the initial conditions of the periodic 
orbit xlv2, the orbits we find have a dynamical behavior similar to the one of the orbit "u" (section 3.1). 
By comparing the orbital behavior of all orbits having a morphology like "u" , with the form of the unstable 
manifolds we conclude that the initial conditions of "u" are very close to the unstable eigendirection of 
the manifold. Both the consequents of "u" and the asymptotic curves wind around the rotational tori 
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Fig. 15. The first 2580 consequents of the orbit "u" (with blue color) sticky to the asymptotic curves of the part A of the 
unstable manifold (with red color) that surround the xl rotational tori, (a) in the 3D subspace (x,x, z). Our point of view 
in spherical coordinates is (0, (ft) = (78°, 84°). (b) in the 3D subspace (x,z,z) with point of view in spherical coordinates 
(0» = (74°,82°) 

associated with xl. The consequents of "u" follow the asymptotic curves and stick in the region dictated 
by the manifolds. For example in Figs. [T5h and [15b the first 2580 consequents of the orbit "u" (depicted 
by blue color) stick close to the asymptotic curves of part A of the unstable manifold (drawn with red 
color). This number of consequents correspond to a time of 124 periods of rotation of our system, a very 
large time interval in galactic dynamics (of the order of several Hubble times). 

For a larger number of consequents, the points depart from part A of the unstable manifold. The relation 
between the consequents of "u" and the unstable manifold is best understood again in the 3D projection 
(x, z, z). As we can see in this projection in Fig. [T6l the first 1000 consequents, after the departure of the 
consequents from part A of the unstable manifold, remain close to the asymptotic curves of part B of 
the unstable manifold. This is the mechanism that creates in phase space the lobes of the ® and the ring 
structure we observe as denser regions in the (z, z) projection in Fig. [5j These denser regions are created 
due to the phenomenon of stickiness around the rotational tori of xlvl and xlvl'. The lobes of the ® 
structure, indicated with white arrows in Fig. [5l are formed by consequents that stay almost on the (red) 
asymptotic curves of Fig. [TBI However, even after the departure of the consequents from the lobes they 
stay close to the xlvl and xlvl' rotational tori, because they are trapped in the immediate neighborhood 
of the eight islands of stability we observe as small black line segments in Fig. [5j In this projection the 
figure resembles what is found in 2D systems close to the last KAM curve (Contopoulos 2002). However, 
we should have in mind that in the 4D space of section of our 3D Hamiltonian system consequents can be 
found on both sides of a torus, moving through the additional dimensions. 

In section 4.2 we said that for the values of perturbation listed in Table 1, the orbits go directly around 
the rotational tori of xlvl and xlvl'. All these orbits have also a sticky behavior at the borders of the 
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Fig. 16. The first 1000 consequents of the orbit "u", given with green "+" symbols after their departure from part A of the 
unstable manifold, as observed in the 3D projection (x, z, z). These consequents stick close to the asymptotic curves of part B 
of the unstable manifold (drawn with red color). Our point of view in spherical coordinates is (6, cj>) — (78°, 79°). 

xlvl and xlvl' 4D stability islands, which we can always see best in the 3D projection (x, z, z) and the 2D 
projection (z,z). In all other projections the sticky regions are discernible as denser regions of the phase 
space embedded in clouds. The (x, z, z) projection has the advantage that in a certain range of viewing 
angles only a few consequents from the cloud are projected in the regions, which look "empty" (actually 
they are occupied by the rotational tori around the stable orbits). 

Our study of the dynamical behavior in the neighborhood of the simple unstable orbits, shows that 
the characteristic double loop structures that have been already found by Magnenat (1982) and Patsis and 
Zachilas (1994) etc., as well as the Q>shaped structures and rings we observe in Fig. [5j are associated with 
stickiness at the borders of rotational tori of stable families that exist in the system. The consequents follow 
the manifolds that wind around the tori. Thus it would be more accurate to associate these structure of 
the phase space mainly with the {7—7-5 transitions, as they are not expected to exist away from such a 
transition point or after a DU —> U transition. 

5. Double Instability 

Double instability in rotating galactic potentials usually appears for large energies. In the xl orbital tree 
(Skokos et al 2002) DU appears in higher order bifurcations, when a stable 3D bifurcation of xl becomes 
first simple unstable and then double unstable, or when a simple unstable bifurcation of xl becomes double 
unstable. Frequently we find DU p.o. at energies for which nearby orbits escape. Thus the selection of the 
appropriate case for investigating the phase space structure is crucial. Here we consider first the dynamics 
in the neighborhood of orbits of the xlv2 family in the region DU (Fig. [1]). However, in order to examine 
the dependence of the dynamical behavior on the energy of the orbit, we study also the phase space close 
double unstable z-axis orbits. These are p.o. along the rotation axis z in rotating triaxial potentials (Heisler 
et al. 1982). 



The periodic orbits of the xlv2 family are DU for large energies (Fig. [I]). At first we study the case for 
Ej = -3.38. The periodic orbit has initial conditions (x , x , z , z ) = (0.088998826,0,0,2.3614474). By 
adding even a tiny perturbation (e.g. Ax = 10~ 6 , 10 -5 orlO -4 ) to its initial conditions, along any direction, 
we find clouds of points in all 3D projections of the 4D space of section These are 4D clouds of points since 
we find scattered points in all 3D projections and mixing of colors when we apply the method of color 
and rotation. Their consequents occupy a volume in the 4D space of section, lying roughly in the region 
(x 1 ,x 2 ) x (xi,x 2 ) x (zi,z 2 ) x (ii,i 2 ) = (-4,4) x (-3,3) x (-2.5,2.5) x (-3,3). 

Then we explore the structure of the phase space in the neighborhood of the double unstable periodic 
orbits of xlv2 in the region DU from —3.41 < Ej < —2.3, for the same fixed value of the perturbation, 



5.1. xlv2 
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Ax = 10~ 4 . We find again 4D clouds of points as in the case of the double unstable periodic orbit at 
Ej = -3.38. 

We study the morphology of the asymptotic curves of the unstable manifold in the neighborhood of the 
double unstable periodic orbits in the region DU and we present this morphology in the neighborhood of 
the double unstable periodic orbit for Ej = —3.38 (Fig. [T7|) . The computation of the asymptotic curves is 
similar to that of the case of the simple unstable periodic orbits (section 4.3). The only difference is that we 
have now two unstable eigenvectors. These are (x±, x±, z\, z\) = (0.0216618,0.163481,-0.737238,0.655198) 
and (x 2 ,x 2 ,,z 2 ,i 2 ) = (-0.831661,0.523481,0.0173103,0.184412). The corresponding eigenvalues are Ai = 
69.08801 and A2 = 1.74981. The initial conditions (x^, X3, Z3, £3), that we use for the computation of the 
asymptotic curves, are given by: 

X3 = X + C\X\ + C 2 X 2 
X3 = X + C\X\ + C 2 X 2 

Z 3 = Zq + C\Z\ + c 2 z 2 

z 3 = z + ciii + c 2 i 2 

We computed the asymptotic curves, in two cases which we named "I" and "II" respectively. In the 
first case, we kept c 2 constant and equal to zero and we used c\ as a parameter. In the second case, we 
kept ci = and we varied c 2 . In case "I" we have taken values c 2 = and 10 -7 < c\ < 2 x 10~ 4 with 
a step 10~ 7 (2000 initial conditions). In the second case "II" we took values c\ = and we varied c 2 in 
the interval -10~ 5 < c 2 < -1.35 x 10 _1 with a step -10~ 5 (13500 initial conditions). In both cases "I" 
and "II" we computed only 3 consequents for every initial conditior0. Here we present the morphology 
of the asymptotic curves in the 3D (x,x,z) projection. We found in our numerical experiments, that the 
morphology of the asymptotic curves in all 3D projections is similar. As it was expected the consequents 
along the asymptotic curves in case "I" (with the large eigenvalue Ai) move away faster from the periodic 
orbit than the consequents along the asymptotic curves in case "II" (with the small eigenvalue A 2 ). 

The consequents of the unstable invariant manifold in case "I" depart from the periodic orbit by 
following the black arrows of Fig. [T71 They form first a loop in the 3D projections of the 4D space of section 
below the periodic orbit. Then the consequents approach again the periodic orbit but with a different color 
(yellow) than the one of the initial part which is violet. Afterwards the asymptotic curve performs a larger 
loop before leaving the box of the frame. In addition we observe a smooth color distribution from violet 
(close to the p.o.) to blue, green, yellow and red and then follows a color oscillation from red to green and 
back again. These oscillations correspond to oscillations of z in the upper part of the color spectrum (right 
part of Fig. \T7\i . This smooth color variation along the curve indicates that these asymptotic curves are 
4-dimensional curves. 

In both cases "I" and "II" the asymptotic curves deviate soon from the eigendirections. In Fig. [18] 
we compare the asymptotic curves "I" (blue color) and "II" (red color) that start along the unstable 
eigendirections that correspond to the large (Ai) and small (A 2 ) eigenvalues respectively. The morphologies 
of these curves are different. The orbits we start integrating close to these curves remain sticky to them 
only for a short time (0.4 periods of rotation of our system) before moving away from them. Soon they 
fill in a chaotic way a large region of the phase space. We have found similar morphologies for all cases of 
stable and unstable manifolds we have examined in the DU region. 

This is a dynamical behavior typical for DU periodic orbits of the xl tree. However, it is not necessarily 
typical for the dynamics in the neighborhood of DU periodic orbits in general. We investigated the orbital 
behavior close to double unstable p.o. in several cases in our system. Below we describe the dynamics close 
to z-axis DU p.o.'s, which differs in many aspects from what we have found in the neighborhood of the 
DU xlv2 p.o.'s we presented up to now. 



1 During the computation of asymptotic curves, we connect gaps that appear along them by considering more initial conditions 
in these gap regions. 



(4) 
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Fig. 17. The asymptotic curves of the unstable invariant manifold in case "I" , for 3 consequents of 2000 initial conditions of 
the initial segment with ci — (a total of 6000 consequents), in the neighborhood of the double unstable periodic orbit for 
Ej — —3.38 in the 4D space of section. The consequents are depicted in the (x,x,z) subspace and are colored according to 
their z value. The location of the p.o. is indicated with a long arrow. The consequents depart from the p.o. in the direction 
pointed by arrows. The point of view in spherical coordinates is (8, <f>) — (45°, 180°). 




Fig. 18. The 3D subspace (x, x, z) of the asymptotic curves of the unstable invariant manifold for case "I" (blue color) 
and case "II" (red color), in the neighborhood of the double unstable periodic orbit for Ej = —3.38 in the 4D space of 
section. The two unstable eigendirections that correspond cases "I" and "II" are depicted with a magenta and a green line 
respectively. The location of the periodic orbit is given with a black square symbol. The point of view in spherical coordinates 
is (9,<P) = (63°, 168°). 



5.2. z-axis 

As already mentioned, the orbits of the z-axis family lie entirely on the rotational, z axis, of the system. 
Thus, we consider now the z = surface as surface of section with z > 0. Earlier studies have shown that 
this family becomes double unstable in slowly rotating triaxal systems (Martinet & de Zeeuw 1988, Patsis 
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Fig. 19. (a) The 2D (a;, A) projection of the orbit we study in the neighborhood of the DU z-axis orbit at Ej = —4.9. We 
consider here 4200 consequents. The first 250 are plotted with blue color. The next 3950 form the red "8"-shaped structure. 
Then they start to diffuse in phase space. Black symbols indicate the location of the periodic orbits in this projection, (b) 
The first 4200 consequents are depicted now with heavy red points, while the next 26000 (light red dots) diffuse and occupy 
all available phase space. 



& Zachilas 1990). We investigated the stability of this family in the potential (2) rotating very slowly with 
Qb = 5 x 10 -6 . For this value of fib the z-axis family has large regions of double instability (DU). In this 
case, as we vary the energy, the z-axis family for large negative values of Ej is stable and for Ej = —6.7 it 
becomes simple unstable. At this transition bifurcates the stable family of p.o. of stable anomalous orbits, 
abbreviated to sao (Heisler et al 1982), and its symmetric family with respect the yz-plane. Then the z-axis 
family, for Ej = —5.836, becomes double unstable and it remains DU for larger values of Ej. 

In the (x,x,y,y) 4D space of section the initial conditions of the z-axis orbits are (0,0,0,0). For 
Ej = —4.9 we present the case of the orbit in the neighborhood of the z-axis p.o., which deviates from it 
in the x direction by Ax = 0.62. We call this orbit "duzl". 

In the 2D (x, x) projection we observe the first 250 consequents being distributed around the p.o. (blue 
points in Fig. [19k). These points do not show clearly any particular structure. The next 3950 consequents 
however, form an "8"-figure. We plot with red color 950 of them in Fig. 119b,. In the same figure the periodic 
orbits existing in this energy are given with black symbols. In the center of the "8" -figure structure the 
x symbol indicates the DU z-axis orbit, while in the central regions of the two lobes of this structure 
the black dots (on the left and the right) indicate the location of sao and its symmetric p.o., which are 
surrounded by projections of rotational tori around them. If we continue to integrate the orbit, we find 
that after more than 4200 intersections the consequents diffuse and occupy all available phase space. Then 
the (x, x) projection is as we see in Fig. [19b . 
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Fig. 20. The time evolution of the finite time LCN of the orbit in Fig. 1191 The arrow in (a) indicates the time at which 
the consequents depart from the "8" -shaped structure that they formed during the first 4200 intersections. After that point 
LCN(t) increases as we can see in (b). 



Table 2. The range 

of perturbations that give orbits whose 
consequents are distributed as in Fig. 1191 



-0.2 < Ax < 0.2 


-1.5 <A±< 1.5 


-0.7 < Ay < 0.7 


-1 < Ay < 1 



On one hand we have the formation of a structure discernible in projections of the phase space. On the 
other hand, the consequents that form the "8"-figure structure of Fig. 119b are not on a smooth 4D surface. 
The application of the color-rotation method shows that we can hardly observe a smooth color variation 
along any direction. These are indications of a chaotic behavior. 

The finite time LCN index points also to the chaotic character of the orbit. For most of the time when 
the consequents in the (x,x) (or the (x,x,y)) projection stay on the "8"-shaped figure, the index varies 
around 5 x 10~ 2 and finally decreases until the time pointed by an arrow in Fig. [20b . This time corresponds 
to 4200 consequents. In Fig. 120b we see how LCN(t) evolves as the consequents start filling all available 
phase space. It increases continuously up to t = 28500 and beyond. 

The range of perturbations for which we find the "8"-figure is registered in Table 2. Outside the range 
given for the x coordinate we find orbits sticky to rotational tori around the sao and its symmetric orbit. 
For absolutely larger perturbations in the x, y or y-direction we encounter three types of orbits, namely (a) 
orbits represented by clouds of points, (b) orbits sticky to tori around the stable periodic orbit xlvl and 
(c) orbits sticky to chains of rotational tori belonging to p.o. of higher multiplicity. This orbital behavior 
is related with the unstable manifold of the z-axis orbit as follows. 

The computation of the manifolds is similar to that of the case of the double unstable periodic orbits 
of xlv2. The unstable eigenvectors are (a?i, x x , yi, yi) = (-0.20427,0.000957614,0.978912,-0.00197504) 
and (x 2 ,x 2 ,y 2 ,y 2 ) = (-0.000880093,0.436429,0.00422376,-0.899728). The corresponding eigenvalues are 
Ai = 5.93858 and A2 = 5.90405. The initial conditions (X3, X3, 2/3, 1/3), which we use for the computation of 
the asymptotic curves, are given by: 

X3 = X + C1X1 + c 2 x 2 
Xz = Xq + ci± 1 + c 2 x 2 

( 5 ) 

2/3 = Vo + cm + c 2 y 2 
m = ya + cil/i + c 2 y 2 

We calculate the manifold by using for c\ the values — 10 -3 < c\ < 10~ 3 with a step 10~ 5 and for 
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Fig. 21. (a) The 4D representation of the unstable asymptotic surface, using the 3D projection (x,x,y) for the spatial 
representation of the consequents and y for the normalized color values. The point of view in spherical coordinates is (8, <f>) — 
(20°, 190°). (b) The (x,x) projection showing the oscillations of the asymptotic surface close to the z-axis DU p.o. located at 
(0,0). 

every value of c\ we consider for C2 the values — 10~ 3 < C2 < 10 -3 again with a step 10 -5 . For every initial 
condition we compute 10 consequents. Then we have totally 4 x 10 5 consequents. This way we computed 
the unstable asymptotic surface we present in Fig. [2TJ We note that now we have an asymptotic surface 
and not an asymptotic curve. This is expected since we have two eigenvalues outside the unit circle (see e.g. 
Arnold 1988, p. 287). The asymptotic surface is easier to be visualized than in the case of the orbit in the 
neighborhood of the double unstable xlv2 orbit, because now the two associated eigenvalues are very close. 
We observe that the unstable asymptotic surface we calculated has a morphology of an "8" -figure in the 
3D projection (x, x, y) of the surface of section. For the 4th dimension we color the points according to the 
normalized i values. We have a complicated but smooth color variation on the surface of the "8" -shaped 
structure (Fig. [2Tb) . In Fig. 121b we give the (x,x) projection of the same figure that shows clearly the 
central oscillations of the surface close to the p.o. located at the center of this "8" -figure structure. 
A similar morphology is found for the stable asymptotic surface. 

The structure we present in Fig. [21] is part of the unstable manifold. If we consider a larger number 
of consequents for —1 < cl,c2 1 and a step 1CP 3 in (5), the asymptotic surface expands considerably 
also in the y-direction and its (x, x) projection is as the red points show in Fig. [22a . We realize that all 
consequents of "duzl" (black color) stick on the unstable asymptotic surface. 

We have to note also that inside the lobes of the "8" -shaped manifold structure, besides the sao 
rotational tori (green elliptical curve in Fig. [22b) we find also sticky orbits, initially forming toroidal 
objects in this region and then diffusing to a larger volume, remaining however sticky to the surface of the 
manifold. The consequents of two such orbits are given with blue and magenta colors in the right lobe of 
the manifold in Fig. [22b . The existence of the asymptotic surface decelerates the diffusion of these orbits 
to larger phase space regions for long times. 

6. Diffusion in phase space 

The study of the orbits in the neighborhood of simple and double unstable periodic orbits gives us an 
information of practical interest namely their rate of diffusion and the possible trapping of their consequents 
in particular regions of the available phase space. This would allow us to separate the orbits that support 
physical structures from those which do not. For this reason we measure the diffusion of the orbits we 
study and we compare the results among themselves. 

From the orbits close to xlv2 p.o. we have consider the following orbital behaviors: 

• The first type of orbits we found close to xlv2 periodic orbits is represented by double loop ribbons or 
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Fig. 22. (a) The (x,x) 2D projection of the unstable asymptotic surface (red color) and the consequents of orbit "duzl" that 
remains for long times sticky to it. (b) The same projection of the unstable asymptotic surface, together with two orbits (blue 
and magenta consequents) that we find by perturbing orbits of the sao rotational tori (projected to a small green elliptical 
curve near the center of the right part of the manifold). 

by closed ribbons in the 3D projections of the 4D space of section. An example is the orbit "u" presented 
in section 4.1. (Figs. \'1\?>\) . The orbits of this orbital type continue beyond the initial ribbons, by forming 
additional dense regions in phase space after long integration times. These dense regions are discernible in 
some projections as "8" -like structures or rings (Fig. [5j). 

The second type of orbital behavior is represented by clouds of points in the 4D surface of section. As an 
example, we consider the orbit with Ax = 10 -4 from the initial conditions of the simple unstable xlv2 for 
Ej = —A and we call it "ul" . 

The third type of orbital behavior is represented again by clouds of points in the 4D surface of section, but 
this time in the neighborhood of a double unstable periodic orbit. As a typical example for the discussion 
to follow, we consider the orbit deviating from xlv2 by Ax = 10 -4 at Ej = —3.38. We call this orbit "dul". 

In addition we include an orbit in the neighborhood of a z-axis DU p.o., i.e. 

The orbit that builds an "8" -shaped figure before diffusing and occupying all available volume of the phase 
space (Fig. [19]) . which we call "duzl". 

In order to estimate the diffusion velocity for the above types of orbital behavior we define a diffusion 
time and a diffusion velocity as follows: In general the volume of phase space occupied by the consequents 
of a chaotic orbit increases with time until they occupy the maximum available volume restricted by the 
surface of the zero velocity. The time that the orbits need to occupy the maximum volume V max , we define 
as "diffusion time" Tj. Then the "diffusion speed" is the ratio 

Vmax 

The quantities V max and Td for the orbits in the neighborhood of a periodic orbit can be easily computed 
in two steps: 

(1) As time increases we compute the mean distance from the periodic orbit s up to this time t of all 
consequents and then we compute a mean volume V(t). In other words we compute a "mean volume" 
of the hypersphere centered on the periodic orbit with radius the mean distance s, i.e the quantity: 
V(t) = ±vrV. 

(2) After time the consequents occupy the maximum volume V max . This is represented by the maximum 
peak in a V(t) diagram and this allows us to compute the diffusion speed. 

In Fig. [23k we plot the volume occupied by the consequents of the orbit "u" (Figs. I2I3D versus time t 
for an interval t < 28000, corresponding to 4200 intersections. We observe that for t < 17300 in our system 
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Fig. 23. The mean volume V(t) that is occupied by the consequents of orbit "u" in the neighborhood of the simple unstable 
periodic orbit for Ej — —4.66. (a) t < 28000, corresponding to 4200 consequents. The arrows indicate the times corresponding 
to 2580 and 3130 intersections, (b) The evolution of V(t) up to t = 7 x 10 5 . 



units (indicated by the left arrow in Fig. [23k corresponding to 2580 intersections) V(t) remains close 
to 0.165 after some initial oscillations. This time interval corresponds to the phase that the consequents 
stick to the asymptotic curves of part A of the unstable manifold forming the double loop ribbon (section 
4.1). Then we observe a second interval between the two arrows in Fig. [23k . where V(t) increases but 
the slope of the curve is smaller than the slope beyond the right arrow. The second arrow corresponds to 
time t = 21800 or 3130 intersections. The dynamical behavior between the two arrows corresponds to the 
time during which the consequents form the lobes of the CD-structure and the "ring" in the (z, z) projection 
(Fig. [5]) sticking on the asymptotic curves of part B of the unstable manifold (section 4.4). Beyond the right 
arrow V(t) increases continuously with a characteristically steeper slope than before. The consequents leave 
all structures observed in the neighborhood of xlv2 (section 4.1). Figure [23b is a continuation of Fig. [23k 
and shows that the peak of V(t) corresponds to V max = 13.433 for = 510000. The speed of diffusion is 
Ud = ^° = 2.63 x 10~ 5 . After the peak of V(t) in Fig. 123b the consequents arrive at distances from the 

periodic orbit that are smaller or equal to the maximum distance. This means that V(t) decreases a little 
and finally, it levels off to a constant value. This happens for all chaotic orbits. 

The evolution of V(t) and the quantity change in the case of the orbit "ul". In Fig. [2"H we see 
that for the curve labeled "1", corresponding to the orbit "ul", we have V ma x = 117.397 and Td = 3180. 
This means that the speed of diffusion is now Ud = 0.0369, i.e. it is 1403 times larger than the velocity 
of diffusion of orbit "u" . The upper, red curve, labeled "2" , describes the evolution of V(t) for the orbit 
"dul" in the neighborhood of the DU p.o. xlv2 at Ej = —3.38. We observe that the curve has a maximum 
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Fig. 24. The evolution of V(t) for the orbit "ul" in the neighborhood of the simple unstable periodic orbit xlv2 for Ej — —4 
(1) and for the orbit "dul" (2), when xlv2 became double unstable for Ej — —3.38. 

V-max = 458 and T<i = 3400. In this case the speed of diffusion is Ud = 0.135 that is 3.65 times and 5133 
times larger than the diffusion speeds of the orbits "ul" and "u" respectively. For "ul" and "dul" V(t) 
increases rather fast to the peak of the diagram, indicating that we do not have any important stickiness 
in these cases. 

Finally for the case we have found stickiness in chaos close to a DU z-axis p.o., i.e. for the orbit "duzl", 
the evolution of the quantity V(t) is given in Fig. [25j Initially V(t) decreases to a value about 1.68 during 
the first 1600 consequents. This period ends at the time indicated by arrow "1" . This happens because the 
majority of the points occupy the right lobe of the "8" structure, as in the (x,x) 2D projection (Fig. [19]) . 
This implies that the majority of the points are included in a volume smaller than the volume occupied 
by the consequents forming the whole "8" -shaped structure. Thus, the mean volume decreases. Then the 
next consequents complete the "8" structure by populating both lobes with similar number of consequents, 
so V(t) increases until V(t) = 3.11969, which corresponds to 3200 intersections (arrow numbered "2"). 
Then we observe a small plateau until the point indicated by the arrow labeled "3" (V(t) = 3.02259). The 
time at this point corresponds to 4200 consequents. This is the number of consequents that stay on the 
eight-figure surface before they depart from it and occupy a larger volume of the phase space. For larger 
times V(t) increases. In Fig. [25b we have a peak at V max = 15.77 for = 1.88 x 10 5 . The diffusion speed 
is Ud = 8.385 x 10" 5 . 

It is evident by comparing Figs. I23ll24l and[25 | that the diffusion of chaotic orbits from the neighborhood 
of simple or double unstable p.o. does not depend on the kind of instability of the periodic orbit. The 
phenomenon of stickiness that decelerates the diffusion appears both to orbits starting close to simple as 
well as to double unstable periodic orbits. When this happens we observe plateaus in the V(t) variation. 

7. Conclusions 

In this paper we studied the structure of the phase space close to simple unstable and double unstable 
periodic orbits in a triaxial 3D system representing a rotating disk galaxy. We have encountered different 
orbital behaviors and different diffusion rates from the neighborhood of the periodic orbit. The phenomenon 
of stickiness seems to be ubiquitous and does not depend on the kind of instability of the periodic orbit 
in the center of the region we study. Our results refer to orbits in the neighborhood of the families xlv2 
(Skokos et al 2002) and z-axis (Martinet & de Zeeuw 1988) and can be summarized as follows: 

• The orbits starting in the neighborhood of a simple unstable periodic orbit form (a) double loop ribbons 
in the 4D space, (b) structures discernible in some projections after long integration times, (c) clouds of 
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points around the U p.o. with mixing of colors in their 4D representations. The two first kinds of orbital 
behavior appear close to the transition U — > S of the xl family at the point of bifurcation of the xlv2 
family and indicate the presence of the phenomenon of stickiness in the dynamics of the system. 
The asymptotic curves of the U xlv2 p.o. close to the U —> S transition of the central family of periodic 
orbits xl are composed of two parts. Part A (calculated for c < in (3)) winds around deformed rotational 
tori belonging to orbits deviating from the xl orbits in the dbz or ±i-direction. Part B (calculated for c > 
in (3)) winds around the rotational tori of the symmetric xlvl and xlvl' p.o. 

In both cases chaotic orbits are guided by the manifolds close to the regions occupied by the rotational 
tori associated with the stable periodic orbits existing in the region. These chaotic orbits are sticky and 
play a very important role, especially in Galactic Dynamics, since the stickiness time is in many cases of 
the order of a Hubble time or even larger. 

Double loops and CD-shaped structures surrounded by rings that are found in several 2D and 3D projections 
of the 4D spaces of section in 3D Hamiltonian systems are formed by orbits close to the asymptotic curves 
of simple unstable periodic orbits that surround the rotational tori of two nearby stable p.o. as in Fig. [5l 
Such configurations appear if the parent family (in our case the planar xl family) has a S — > U transition 
and a nearby U —> S transition. At the first transition bifurcate two symmetric stable orbits (like xlvl and 
xlvl') and at the second transition we have the bifurcation of two simple unstable families (in our case the 
families xlv2 and xlv2'). Near the U — > S transition of the parent family, the already bifurcated stable 
orbits are not far from the unstable orbit under study, thus the asymptotic curves of the latter families 
can surround their rotational tori. Away from the transition points of the varying parameter (in our case 
the Jacobi constant) the shape of the asymptotic curves become complicated and the consequents of the 
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nearby orbits form clouds in the 4D spaces of section. 

Double instability in 3D rotating galactic disks appears for larger energies than simple instability in the 
families of the xl tree. The usual case is after a U — > DU transition and in general there are no rotational 
tori associated with stable families in the neighborhood of the DU p.o. Typically the consequents of the 
orbits form clouds of points in the 4D spaces of section. 

In the case of a double unstable periodic orbit there are two different pairs of eigendirections, starting 
at the periodic orbit. The curves of each pair start in opposite directions. Between the eigendirections an 
"eigensurface" is formed. However the shape of this surface is complicated and it is difficult to be visualized 
when the associated eigenvalues are not close to each other. 

In the case of z-axis DU periodic orbits in slow pattern speed models, at low energies, we find that the 
eigenvalues corresponding to the unstable eigenvectors are close to each other. Because of that the speeds 
of deviation along the two eigendirections are similar. This allows us to visualize the unstable eigensurface 
corresponding to the DU periodic orbit. We find sticky chaotic orbits whose consequents stay close to the 
unstable manifold for very long times (of the order of a Hubble time or more) . We find these sticky orbits 
either by perturbing the initial conditions of the z-axis DU p.o., or by perturbing orbits on rotational tori 
associated with the stable anomalous orbit (sao) and its symmetric family. 

Among the orbits we studied we found those close to the double unstable orbits of the xlv2 family having 
the largest diffusion speed. The sticky chaotic orbits close to the bifurcating point of the simple unstable 
xlv2 orbit and close to the double unstable z-axis orbit we have examined, have comparable diffusion 
speeds. They are much slower than the diffusion speeds of the orbits in the neighborhood of xlv2 simple 
unstable periodic orbits away from its bifurcation point, or of the double unstable orbits of the same family 
that have eigenvalues of very different size. 

These results will be used in a forthcoming paper to assess the role of chaotic orbits in shaping the 
vertical profiles of disk galaxies. 
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